rng('default');

% simulate to get unconditional joint kappa/lastIssuance density without
% disasters (mimicking the post-war sample)
Binit           = 10000;  % number of initial periods to discart
Ninit           = 100;    % number of firms
Tinit           = 120000;
[~,states]      = sim_markov_chain_fast(P_nodis,Tinit+Binit,1,1,1:nS-1);
g_sim           = mu_X(states) + sig_XS(states).*randn(Tinit+Binit,1) + sig_Xid(states).*randn(Tinit+Binit,Ninit); % [Tsim,Nsim]
k_def_sim       = k_def(states); % [Tsim,1]
k_iss_sim       = k_iss(states); % [Tsim,1]
k_bar_sim       = k_bar(states); % [Tsim,1]
k_sim           = NaN(size(g_sim));
k_QML           = NaN(size(g_sim));
last_issuance   = ones(size(g_sim)); % times initial state of the Markov chain
k_sim(1,:)      = k_bar(states(1));
k_QML(1,:)      = k_bar(states(1));
for t = 1:Tinit+Binit-1 % maturity
    k_lag                = k_sim(t,:);
    G_lag                = exp(g_sim(t,:));
    default              = k_lag <= k_def_sim(t);
    issuance             = k_lag >= k_iss_sim(t);
    k_sim(t+1,:)         = k_lag.*G_lag;
    k_sim(t+1,default)   = k_bar_sim(t).*G_lag(default);
    k_sim(t+1,issuance)  = k_bar_sim(t).*G_lag(issuance);
    k_QML(t+1,:)         = k_QML(t,:);
    k_QML(t+1,issuance)  = k_bar_sim(t);
    last_issuance(t:end,issuance)  = ones(Tinit+Binit-t+1,sum(issuance)).*states(t); % most recent issuance state
end
k_pdf     = k_sim(Binit+1:end,:);
k_QML_pdf = k_QML(Binit+1:end,:);
LI_pdf    = last_issuance(Binit+1:end,:);
k_pdf     = k_pdf(:);
k_QML_pdf = k_QML_pdf(:);
LI_pdf    = LI_pdf(:);
n_pdf     = length(k_pdf);
clear Binit Ninit Tinit states g_sim k_def_sim k_iss_sim k_bar_sim last_issuance k_sim t k_QML



% generate 'disaster' paths
NOBS        = 100*12; % total # of periods to simulate
NOBS_before = 2;  % years before event
NOBS_after  = 10; % years after event
Tsim        = (NOBS_before+NOBS_after)*12; % # of periods saved
Nsim        = 100; % # of firms
Rsim        = 1000; % # of reps
Psim        = 100; % # of paths per rep

% define disaster
Dsize  = -0.1658; % simple growth rate
Dmin   = Dsize-0.0025; % range of cummulative drop
Dmax   = Dsize+0.0025;
D_isD  = -0.05; % no other drop below this value in event window (so that the path contains a single 'disaster')

MOM1   = NaN(Tsim,Rsim);
MOM2   = NaN(Tsim,Rsim);
MOM3   = NaN(Tsim,Rsim);
MOM4   = NaN(Tsim,Rsim);
MOM5   = NaN(Tsim,Rsim);
MOM6   = NaN(Tsim,Rsim);
MOM7   = NaN(Tsim,Rsim);
MOM8   = NaN(Tsim,Rsim);
MOM9   = NaN(Tsim,Rsim);
MOM10  = NaN(Tsim,Rsim);
MOM11  = NaN(Tsim,Rsim);
state1 = NaN(Tsim,Rsim);
state2 = NaN(Tsim,Rsim);
for rr=1:Rsim
    
    % simulate aggregate paths
    delC_paths   = NaN(Tsim,Psim);
    epsE_paths   = NaN(Tsim,Psim);
    state_paths  = NaN(Tsim,Psim);
    i            = 0; % # of crisis paths found
    while 1
        
        % simulate consumption and systematic earnings
        s0          = sum(rand>cumsum(stat(1:nS-1)))+1;
        [~,states]  = sim_markov_chain_fast(P,NOBS,1,s0,1:nS);
        epsCE       = randn(NOBS,2)*chol([1,rho_CE;rho_CE,1]); % correlated shocks for consumption and systematic earnings
        epsC_sim    = epsCE(:,1);
        epsE_sim    = epsCE(:,2);
        delC_sim    = mu_C(states) + sig_C(states).*epsC_sim;
        C_a         = sum(reshape(exp(cumsum(delC_sim)),[12,NOBS/12]))'; % time-aggregate consumption
        delC_a      = [NaN;log(C_a(2:end)./C_a(1:end-1))]; % time aggregate consumption
        
        
        disaster_D  = false(size(delC_a)); % indicator: =1 in first year of disaster
        disaster_all= false(size(delC_a)); % indicator: =1 in first year of disaster with minimum drop of 'D_isD' %
        cum_loss    = 0;
        cum_periods = 0;
        for jj=1:length(delC_a) % search for subsequent years with negative growth
            if delC_a(jj)<0
                cum_loss    = cum_loss + delC_a(jj);
                cum_periods = cum_periods + 1;
            else
                if cum_loss>log(1+Dmin) && cum_loss<log(1+Dmax) % drop of the desired magnitude
                    disaster_D(jj-cum_periods) = true;
                end
                if cum_loss<log(1+D_isD) % some large drop
                    disaster_all(jj-cum_periods) = true;
                end
                cum_loss = 0;
                cum_periods = 0;
            end
        end
        
        % check whether a suitable disaster can be found
        ix = find(disaster_D); % index of year in which the crisis began
        if ~isempty(ix) % at least 1 disaster of desired magnitude ocurred in this sample
            if ix(1)>(1+NOBS_before) && ix(1)<(length(delC_a)-NOBS_after) % sufficiently many observations before and after
                if sum(disaster_all(ix(1)-NOBS_before:ix(1)+NOBS_after))==1 % only one disaster (of at least 'D_isD' %) occured in the event window
                    if ~any(states( (ix(1)-1-NOBS_before)*12+1:(ix(1)+1-NOBS_before)*12 )==3) && states( (ix(1)+1-NOBS_before)*12+1 )==3 % state 3 begins exactly at t=0
                        i                 = i + 1;
                        window            = (ix(1)-1)*12+1-NOBS_before*12:(ix(1)-1)*12+NOBS_after*12;
                        epsE_paths(:,i)   = epsE_sim(window);
                        delC_paths(:,i)   = delC_sim(window);
                        state_paths(:,i)  = states(window);
                    end
                end
            end
        end
        
        if i==Psim
            break
        end
    end
    
    
    % simulate kappa panels
    g_sim                = mu_X(state_paths) + sig_XS(state_paths).*epsE_paths + sig_Xid(state_paths).*randn(Tsim,Psim,Nsim);
    k_def_sim            = k_def(state_paths); % [Tsim,Psim] state-dependent default threshold
    k_iss_sim            = k_iss(state_paths); % [Tsim,Psim] state-dependent issuance threshold
    k_bar_sim            = k_bar(state_paths); % [Tsim,Psim] state-dependent optimal kappa
    init                 = round(rand(Psim,Nsim)*n_pdf+0.5);
    default              = false(Tsim,Psim,Nsim);
    k_sim                = NaN(Tsim,Psim,Nsim);
    k_sim(1,:,:)         = k_pdf(init);
    k_QML                = NaN(Tsim,Psim,Nsim);
    k_QML(1,:,:)         = k_QML_pdf(init);
    last_issuance        = NaN(Tsim,Psim,Nsim);
    last_issuance(1,:,:) = LI_pdf(init);
    
    for t = 1:Tsim-1
        
        kt    = squeeze(k_sim(t,:,:)); % [Psim,Nsim]
        qmlt  = squeeze(k_QML(t,:,:)); % [Psim,Nsim]
        kbt   = k_bar_sim(t,:)'*ones(1,Nsim); % [Psim,Nsim]
        gt    = squeeze(exp(g_sim(t,:,:))); % [Psim,Nsim]
        kt1   = kt.*gt;
        qmlt1 = qmlt;
        kbt1  = kbt.*gt;
        
        default_t        = kt <= k_def_sim(t,:)';
        kt1(default_t)   = kbt1(default_t);
        qmlt(default_t)  = kbt(default_t);
        default(t,:,:)   = default_t;
        
        issuance_t       = kt >= k_iss_sim(t,:)';
        kt1(issuance_t)  = kbt1(issuance_t);
        qmlt(issuance_t) = kbt(issuance_t);
        
        k_sim(t+1,:,:)   = kt1;
        k_QML(t+1,:,:)   = qmlt;
        
        st                     = state_paths(t,:)'*ones(1,Nsim);
        lit                    = squeeze(last_issuance(t,:,:));
        lit(issuance_t)        = st(issuance_t);
        lit(default_t)         = st(default_t);
        last_issuance(t+1,:,:) = lit;
        
    end
    
    
    % quasi-market-leverage
    Fit                = griddedInterpolant({1:nS,k},lamDx);
    DX_sim             = Fit(repmat(state_paths,[1,1,Nsim]),k_QML);
    Fit                = griddedInterpolant({1:nS,k},lamSx);
    SX_sim             = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    QMLEV_sim          = DX_sim./( DX_sim + SX_sim .* k_sim./k_QML);
    QMLEV_sim(default) = NaN;
    
    
    % realized loss rate
    Fit     = griddedInterpolant({1:nS,k},lamD);
    Ds_sim  = Fit(last_issuance,k_QML);
    recov   = lamA.*(1-omega);
    loss    = 1-recov(repmat(state_paths,[1,1,Nsim]))./Ds_sim .* (k_sim./k_QML);
    loss(~default) = NaN;
    
    
    % other firm policies and prices
    Fit         = griddedInterpolant({1:nS,k},EP_R);
    E_R_sim     = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    Fit         = griddedInterpolant({1:nS,k},IV_ATM);
    IV_sim      = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    Fit         = griddedInterpolant({1:nS,k,1:nS},CDS(:,:,:,60));
    CDS60_sim   = Fit(repmat(state_paths,[1,1,Nsim]),k_sim,last_issuance);
    Rf1         = Rf(:,1);
    R_sim       = E_R_sim - Rf1(repmat(state_paths,[1,1,Nsim]));
    Fit         = griddedInterpolant({1:nS,k},PD_q(:,:,60)-PD_p(:,:,60));
    CDS60_spread_sim = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    Fit         = griddedInterpolant({1:nS,k},PD_p(:,:,60));
    PD60_sim    = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    Fit         = griddedInterpolant({1:nS,k},LGD_p(:,:,60));
    LGD60_sim   = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    Fit         = griddedInterpolant({1:nS,k},LGD_q(:,:,60));
    LGD60q_sim  = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    
    
    % stock price index
    default_last_period = false(size(default));
    default_last_period(2:end,:,:) = default(1:end-1,:,:);
    Rex_sim                = NaN(size(R_sim));
    Rex_sim(2:end,:,:)     = exp(log(SX_sim(2:end,:,:))  - log(SX_sim(1:end-1,:,:)) + g_sim(1:end-1,:,:)) - 1;
    Rex_sim(default)       = -1;
    Rex_sim(default_last_period) = NaN;
    R_MKT_sim              = trimmean(Rex_sim,2,3); % expected returns are measured contemporaneously to mktcaps
    stock_price            = NaN(size(R_MKT_sim));
    stock_price(2:end,:,:) = cumprod(1+R_MKT_sim(2:end,:));
    stock_price            = stock_price./(stock_price(NOBS_before*12,:)); % normalize to one in the period where the disaster starts
    
    % save
    MOM1(:,rr)   = 100 *mean(delC_paths,2); %
    MOM2(:,rr)   = 100 *mean(nanmean(CDS60_sim ,3),2); %
    MOM3(:,rr)   = 100 *mean(nanmean(QMLEV_sim ,3),2); %
    MOM4(:,rr)   = 1200*mean(nanmean(default   ,3),2); %
    MOM5(:,rr)   = 100 *mean(nanmean(IV_sim    ,3),2); %
    MOM6(:,rr)   = 100*nanmean(nanmean(loss    ,3),2); %
    MOM7(:,rr)   = mean(nanmean(stock_price ,3),2); %
    MOM8(:,rr)   = 100 *mean(nanmean(PD60_sim ,3),2); %
    MOM9(:,rr)   = 100 *mean(nanmean(CDS60_spread_sim ,3),2); %
    MOM10(:,rr)  = 100 *mean(nanmean(LGD60_sim ,3),2); %
    MOM11(:,rr)  = 100 *mean(nanmean(LGD60q_sim ,3),2); %
    state1(:,rr) = mean(nanmean(repmat(state_paths==1,[1,1,Nsim]) ,3),2);
    state2(:,rr) = mean(nanmean(repmat(state_paths==2,[1,1,Nsim]) ,3),2);
    
end



